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SPIRAL PHASE DEMODULATOR FOR TWO-DIMENSIONAL 

FRINGE PATTERNS 



^ Field of the Invention ^ ~ 

5 The present invention relates to two-dimensional code patterns and more 

particularly to the automatic demodulation of fringe patterns. 



Background 

Methods exist for demodulating fringe patterns. However, conventional methods 
10 of demodulating closed curve fringe patterns, including circular and elliptical fringe 
patterns, inevitably produce ambiguities. Ambiguities are then typically resolved by 
including additional constraints. These additional constraints may include restrictions on 
the smoothness of curves within the fringe patterns. 

In addition to ambiguities other defects and errors may also be produced by the 
15 preliminary demodulation process. An example of this can be found in the classic Fourier 
(Hilbert) transform method (FTM). When the FTM is applied to a whole image, artefacts 
are produced, the artefacts being related to discontinuities introduced with the half-plane 
filter used in Fourier space. The artefacts result in errors in phase estimation, and 
specifically in regions where an angle of the fringe pattern is close to perpendicular with 
20 the half-plane discontinuity line. 

A method adopted to overcome this problem is to choose the half-plane filter in 
Fourier space such that the discontinuity line does not cut through the signal in Fourier 
space. 

However, the frequency components derived from circular or elliptical fringe 
25 patterns are also "circular" and no discontinuity line could be defined which does not cut 
through the frequency signal. 

Other methods based on local (small kernel) demodulation may avoid the errors 
of the FTM, but typically fail to disambiguate correctly when there are strong 
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perturbations in the underlying fringe pattern. Features such as bifurcations and fringe 
endings represent typically strong perturbations. 

Disclosure of the Invention 

5 It is an object of the present invention to substantially overcome, or at least 

ameliorate, one or more disadvantages of existing arrangements. 

According to a first aspect of the invention, there is provided a method of 
demodulating a real two-dimensional pattern, the method comprising the steps of: 

estimating a quadrature two-dimensional pattern from said real two-dimensional 
10 pattern using a two-dimensional spiral phase filter; and 

creating a demodulated image by combining said real two-dimensional pattern 
and said estimated quadrature two-dimensional pattern. 

According to another aspect of the invention, there is provided an apparatus for 
implementing the aforementioned method. 
15 According to another aspect of the invention there is provided an apparatus for 

calculating a quadrature conjugate of a real two-dimensional code pattern, the apparatus 
comprising: 

a first spatial light modulator for modulating coherent light to produce said real 
two-dimensional pattern; 
20 a first lens for Fourier transformation of said real two-dimensional pattern to 

produce a first Fourier transformed image; 

a second spatial light modulator or spiral phase plate for phase modulating said 
first Fourier transformed image to produce a phase modulated image; 

a second lens for Fourier transformation of said phase modulated image to 
25 produce a second Fourier transformed image; 

a third spatial light modulator or spiral phase plate for phase modulating said 
second Fourier transformed image to produce said conjugate of said real two-dimensional 
pattern; and 
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an image sensor for capturing said conjugate of said real two-dimensional 

pattern. 

Brief Description of the Drawings 

A number of 'preferred embodiments of the present invention will now be 
5 described with reference to the drawings and appendix, in which: 

Fig. 1 is an illustration of the local structure of a fringe pattern; 

Fig. 2 is an illustration of the Fourier transform of the local fringe pattern of 

Fig. 1; 

Fig. 3 A is a flow diagram of a method of phase demodulation; 
10 Fig. 3B is a flow diagram of an alternative method of phase demodulation; 

Fig. 4 is an example of a fringe pattern image upon which the preferred 
embodiment of the invention can be applied; 

Figs. 5A and 5B are the real and imaginary parts of an intermediate result in 
applying the method of Fig. 3B to the image of Fig. 4; 
15 Fig. 6 is an estimate of fringe orientation of the fringe pattern image of Fig. 4; 

Figs. 7 A and 7B are the magnitude and phase components of a demodulated 
fringe pattern when the estimate of fringe orientation of the fringe pattern image as shown 
in Fig. 6 is used; 

Fig. 8 is an alternative estimate of fringe orientation of the fringe pattern image 
20 of Fig. 4; 

Figs. 9 A and 9B are the magnitude and phase components of a demodulated 
fringe pattern when the estimate of fringe orientation of the fringe pattern image as shown 
in Fig. 8 is used; 

Fig. 10 is a schematic block diagram of an optical spiral phase demodulator upon 
25 which the preferred embodiment of the present invention can be practiced; 

Fig. 1 1 is a schematic block diagram of a general purpose computer upon which 
an alternative embodiment of the present invention can be practised; and 

Appendix 1 is a paper by the inventors entitled "Generalised Fringe Pattern 
Analysis using a Spiral Phase Transformation". 
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Detailed Description including Best Mode 

Consider a two-dimensional fringe pattern with intensity of the fringe pattern as 

5 follows: 



f(x,y) = a(x,y) + b(x t yXcos[27i(uo{x-x 0 }+v 0 {y-yo})+y/\ (1) 

Position coordinates (jc, y) may be continuous for an analog pattern or discrete 
10 for digital patterns. A background level is denoted by a(x f y) while an amplitude 
modulation term is denoted by b(x f y), both assumed to be slowly varying functions. 
Phase term H/(x 0 ,y 0 ) presents the local offset of the overall phase term in the square 

brackets [] of equation (1). 

Referring to Fig. 1, it is assumed that the frequency components u 0 and v c are 
15 slowly varying so that they are effectively constant over the small region of interest, Q, 
around (jc 0 , yo), making the carrier, a linear function of x and y. 

A noise term is omitted from equation (1), but may be present in real fringe 

patterns due to the occurrence of blurring, non linearities, quantisation errors, smudging, 

scratches, cuts, dust, etc. 
20 The intensity pattern / (jc, y) is (generally) a rapidly varying function of (x, y) . 

The Fourier transform (FT) of the above fringe pattern region Q is: 



F ni u > v )= \\f{x,y)^-2m(ux^vy)}lxdy 
=A o {u,v)^B n {u-u 0 ,v-v 0 )extiiv}+B n {u+u 0 , V +v^^ 



(2) 



25 wherein / = V-T 

From equation (2) it is clear that the Fourier transform of the fringe pattern has 3 
lobe centres, the first being at the origin and another two at (w 0 ,v 0 ) and (-w 0 ,-vo). The last 
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two lobe centres are mirror images of each other with regards to the origin. The effect of 
windowing the fringe pattern over a small region Q is to blur the otherwise sharp lobes. 
The central lobe representing the zero frequency components, can be removed, and 
will therefore be disregarded in the following analysis, '■"Methods' that may be used for 
removal of the zero frequency components include using a suitably matched high pass 
filter. The remaining spectral lobes are illustrated in Fig. 2. 

The local spatial frequency components uq and vo of the fringe pattern satisfies 
the following equations: 



10 



"o + v o = °"o 

tan(A,)=i 



o J 



(3) 



wherein /5b is an angle of the normal to the fringes with the x-axis. 

The conventional Fourier transform method (FTM) of fringe analysis is based on 
the assumption that the angle /5b is confined to a small range for the entire fringe pattern, 

15 typically much less than n radians. However, when considering fringe patterns with 
circular, elliptical carrier patterns or other closed curves, the angle po covers the entire 
range from 0 to 2iz radians. The FTM would therefore be deficient in solving fringe 
patterns created with closed curve carriers. The reason for this failure is related to the 
inability of the FTM to isolate the two (mirrored) lobes of the FT. As noted above, the 

20 frequency components derived from closed curve fringe patterns are also "closed curves". 
Therefore, referring again to equation (2), both the second and third terms will have 
"closed curve" spectral components, causing a total lobe overlap. 

The preferred embodiment of the invention removes directional discontinuities in 
the Fourier processing of the fringe pattern by smoothly encoding the fringe direction by 

25 utilizing a filter with odd (180° rotational) symmetry, but (unlike the FTM) a smooth 
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transition between adjacent frequencies. Also, the filter used in the preferred embodiment 
does not alter the strength of different frequency components, and is therefore a unit 
magnitude or phase only filter. 

A (complex) Fourier representation P(u,v) of the filter preferably has the 
5 following anti-symmetry qualities: 

P(u, v)=exp[z Y{u, v)]=-/>(-w,-v)=-exp[/( v)+re)] (4) 

The preferred parameter chosen has a phase which is set equal to the polar angle: 

10 

*F(«.v) = <f> (5) 



where 



15 



u-q cos(^) 

v = q sin(^) 

2 2 2 
u + v =q 



(6) 



The resulting filter is termed a "spiral phase filter" because the filter phase 
parameter ^increases linearly as the polar angle (j> increases. The magnitude of the filter 
is constant (so the filter is a pure phase filter) and the phase has the fon^ of a spiral or 
20 helical surface with just one turn. An equivalent form of the filter is simply: 

P(M , V) = (^) = ^ (7) 
<7 

The Fourier Transform of the fringe pattern after the central lobe has been 
25 removed is as follows: 

^( U) v)=(5 0 ( M -u 0 ,v-vJexr{iH+^ n ("+"o,v+vJexr{-/^])e-* +%) (8) 
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Applying the spiral phase filter P(u,v) to equation (8) changes the lobe symmetry to gives: 



G(u f v) = P(u,v)^(u,v) 

= (B n (w - u 0 , v - v 0 ) exp[z ^] exp[zf] - 5 n (w + u 0 , v + v 0 ) exp[-z ^] exp[z0])e~ 2 " ( ^ +v:Vo) 

5 

Equation (9) can be inverse transformed, as is shown in equation (10), by 
implicitly assuming that the polar angle ^(w,v) « ± po remains constant over each lobe. 
This is equivalent to assuming that the lobe subtends a small polar angle, which is a 
reasonable assumption for typical fringe patterns. 

10 

g(x,y) = £ £ G(w, v) exp[27u(ux + vy)]dudv ^ 
= -i exp[zp 0 ] sin[27i(w 0 {x - x 0 } + v 0 {y - y 0 } ) + v|/] 



The desired sine component of the local fringe pattern g(x,y) can then be 
extracted by multiplying g(x,y) with an orientational phase component, exp [-//%] as 
15 follows: 

4 

g(x,y).exp[-ip 0 ] =-isin[27V (u 0 {x - x 0 } + v 0 {y - y 0 }) + ^7 (11) 



Conceptually the process is repeated for all fringe regions with their 
20 corresponding fringe orientation angles J3. In practice the procedure can be applied to all 
regions simultaneously. It is noted that the original (conceptual) windowing is not 
actually required. 

The aforementioned process depends upon the functions b(x.y), p(x,y), y/(x,y), 
Uo(*,y), v 0 (x,y), and <J 0 (x,y) having suitably slow variation. However, it is found in 
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practice that the preferred embodiment can be successfully implemented on any fringe 
patterns with suitably smooth variation in orientation and spacing. For fringe patterns 
with circular fringes and constant or quadratic (chirp) fringe spacing the preferred 
embodiment only breaks down near the very center of the fringe pattern due to the 
5 assumption of small polar angle subtense failing near the centre. 

Equation (11) represents an estimation of the imaginary part of the analytic 
function (or image) related to the original real fringe pattern. This component is also 
called the quadrature component or analytic conjugate. By forming a new (complex) 
output image comprising the original real image and the estimated imaginary image, the 
10 phase and magnitude of the resulting complex fringe representation can be obtained as 
follows: 

if ' a)+{g exp{-z/%» « b.exp{-i[2m(u 0 {x-XQ}+vo{y-y 0 })+i//]} (12) 

15 A flowchart showing the phase demodulation process is illustrated in Fig. 3 A. 

The fringe pattern image f(x,y) 9 captured in step 10, typically includes only a real part. 

The fringe pattern image f(x,y) is transformed in step 20 by applying a 2-dimensional 

FFT. Because the fringe pattern image f(x,y) is real, the real and imaginary parts of the 

transform F(u, v) are symmetric and antisymmetric respectively. 
20 The DC component (lobe) of the transform F(u,v) is removed in step 30. In step 

40, the spiral phase filter P(u,v) is applied to the high passed filtered transform F(u,v). 

The inverse Fourier transform is calculated in step 50. This may be performed by again 

applying a 2-dimensional FFT. 

In step 60, an orientational phase filter is applied to the resulting image of step 
25 50. In step 70 the real part is discarded, leaving only the imaginary part, which is then 

added to the input fringe pattern image f(x,y), in step 80 to produce an output image in 

step 90. 

A flowchart showing an alternative embodiment of the phase demodulation 
process is illustrated in Fig. 3B and will be explained with reference to an example fringe 
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pattern. An example of such a two-dimensional fringe pattern is illustrated in Fig. 4 
wherein only the real part of the pattern is known. 

The fringe pattern image f(x,y) is captured in step 210. In the alternative 
embodiment of the phase demodulation process, steps 20 to 50 of the process in Fig. 3 A 
5 are replaced by a single convolution step 220. The fringe pattern image f(x,y) is 
preferably convolved with spiral phase function p(x,y), wherein p(x,y) is defined as 
follows: 



-yfu 

10 



P( x >y) ~ I I P(u,v)exp[2m(ux + vy)]dudv 

; r exp[2m(ux + vy)]dudv 

CO J-OO f_ 2 I y2 



= '-^P- = T^-expOe) (14) 

2nr 2nr 



wherein 

x = rcos^| (15) 

y = r sin 

15 The result after such a convolution on the fringe pattern image f(x,y) of Fig. 4 is 

illustrated in Figs. 5 A and 5B showing the real and imaginary parts of f(x,y)*p(x,y). 

Referring again to Fig. 3B, in step 230 the orientational phase filter, exp[-i(J] is 
applied. However, an estimation of the fringe orientation |3 is needed. 

A method of estimating the fringe orientation p is to compute the 2-D energy 

20 operator which gives uniform estimates over the fringe pattern image f(x,y). Typical 
gradient based methods fail to give reliable orientation estimates at the peaks and valleys 
of the fringes. Such a method is described in the publication "A multi-dimensional 
energy operator for image processing", SPBE Conference on Visual Communications and 
Image Processing, Boston, MA, [1992], 177-186, by P. Maragos, A.C. Bouik, and J. F. 
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Quartieri, the contents of which are incorporated herein by cross reference. The 2-D 
energy operator is applied to the DC removed fringe pattern image separately in the x and 
y directions as follows: 

Kx,y)=/tx,y)-a(x,y) (15) 

where h(x,y) is the intensity of the fringe pattern with the DC removed. E x {} and E y {) are 
the x and y energy operators respectively and are defined as folows: 



10 E x {Kx,y)} = \^\ -A0«(2^ o fe) 2 (16) 



E y {h(x,y)} = \^-\ -h^~{2nv 0 b) 2 (17) 



dy) dy 



The 2-D energy operator can be defined as follows: 

15 

E x {h(x,y)} + iE y {h(x,y)}** {27m 0 bf + i{2nv 0 bf = (2nq 0 b) 2 exp(2ij3 0 ) (18) 

An overall ambiguity results, because a fringe pattern oriented at 0° cannot be 
distinguished from a fringe pattern oriented at 180° unless the overall fringe environment 

20 is considered. The ambiguity typically appears as an estimate of twice the orientation 
phase exp[2ip]. However, step 230 requires a exp[-ip] phase factor. A simple square 
rooting gives ±exp[-ip] ? again causing an ambiguity. The variation is with respect to a 
fixed fringe orientation. In Fig. 6 the estimate of the fringe orientation is shown. As the 
fringes pattern crosses a threshold orientation at p=±90°, the ambiguity "flips". This can 

25 be seen in the phase component of a demodulated phase pattern as is illustrated in Figs 
7A and 7B wherein the magnitude and phase components of the result of step 230 is 
shown. Such an ambiguity is tolerable in certain applications. Note that this ambiguity is 
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a sharp "flip" of sign, quite unlike the (error inducing) directional artefacts of the FTM 
demodulated phase. 

An alternative method of estimating the fringe orientation p is to use an uniform 
orientation estimator as in the previous method, but to utilise local environmental 
5 constraints (such as topology of smoothly varying fringes) to resolve or unravel the 
orientational ambiguity. Fig. 8 illustrates a smooth estimate of the fringe orientation 
exp(i/?) derived from exp(2z/?). The orientation phase estimate has only one 

discontinuity line at positions where the phase changes from 360° to 0°. A number of 

different methods are suitable for resolving the orientational ambiguity. The magnitude 
10 and phase components of the result of step 230 using the estimation of the fringe 

orientation (3 of Fig. 8 is illustrated in Figs 9 A and 9B. 

It is noted that an estimate of the orientation can be obtained from the output of 

step 50 in Fig. 3A or alternatively from the output of step 220 in Fig. 3B. However, 

additional unravelling will be required to obtain a useable orientation estimate. 
15 In step 240 the real part is discarded, leaving only the imaginary part, which is 

then added to the input fringe pattern image f(x,y), in step 250 to produce an output image 

in step 260. 

The method of Fig. 3A can be practiced using an optical spiral phase 
demodulator 300, such as that shown in Fig. 10. The optical spiral phase demodulator 

20 300 is illuminated by coherent light 301. A first spatial light modulator (SLM) 302 
modulates the coherent light 301 to produce the input image, which is typically a fringe 
pattern f(x,y). The fringe pattern f(x,y) is Fourier transformed by lens 303. A second 
SLM 304 or a spiral phase plate performs phase modulation in the image projected on it. 
The phase modulated image is then Fourier transformed by lens 305 and is projected onto 

25 a final SLM 306 or spiral phase plate. The final spiral phase plate 306 retards the 
coherent optical beam in a manner such that the plate introduces a spiral phase 
multiplication. An image sensor 307 is in close proximity to the final SLM 306 to detect 
an output image. The output image will be in quadrature (or analytic conjugate) to the 
input image. In other words, if a cosine pattern is input, a sine pattern is output. 
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In the special case where the input image is a circularly symmetric fringe pattern, 
then the SLM 306 has the simple form of a spiral phase with opposite sign to SLM 304. 

The method of Fig. 3 A or 3B may alternatively be practiced using a conventional 
general-purpose computer system 100, such as that shown in Fig. 11 wherein the 
5 processes of Fig. 3A or 3B may be implemented as software, such as an application 
program executing within the computer system 100. In particular, the steps of method of 
phase demodulation are effected by instructions in the software that are carried out by the 
computer. The software may be stored in a computer readable medium, including the 
storage devices described below. The software is loaded into the computer from the 

10 computer readable medium, and then executed by the computer. A computer readable 
medium having such software or computer program recorded on it is a computer program 
product. The use of the computer program product in the computer preferably effects an 
advantageous apparatus for demodulating two-dimensional fringe patterns in accordance 
with the embodiments of the invention. 

15 The computer system 100 comprises a computer module 102, input devices such 

as a keyboard 110 and mouse 112, output devices including a printer 108 and a display 
device 104. 

The computer module 102 typically includes at least one processor unit 114, a 
memory unit 118, for example formed from semiconductor random access memory 

20 (RAM) and read only memory (ROM), input/output (I/O) interfaces including a video 
interface 122, and an I/O interface 116 for the keyboard 110 and mouse 112. A storage 
device 124 is provided and typically includes a hard disk drive 126 and a floppy disk 
drive 128. A magnetic tape drive (not illustrated) may also be used. A CD-ROM 
drive 120 is typically provided as a non-volatile source of data. The components 114 

25 to 128 of the computer module 102, typically communicate via an interconnected bus 130 
and in a manner which results in a conventional mode of operation of the computer 
system 100 known to those in the relevant art. Examples of computers on which the 
embodiments can be practised include IBM-PC's and compatibles, Sun Sparcstations or 
alike computer systems evolved therefrom. 
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Typically, the application program of the preferred embodiment is resident on 
the hard disk drive 126 and read and controlled in its execution by the processor 114. 
Intermediate storage of the program may be accomplished using the semiconductor 
memory 118, possibly in concert with the hard disk drive 126. Tn some instances, the \ 
5 application program may be supplied to the user encoded on a CD-ROM or floppy disk 
and read via the corresponding drive 120 or 128, or alternatively may be read by the user 
from a network via a modem device (not illustrated). Still further, the software can also 
be loaded into the computer system 100 from other computer readable medium including 
magnetic tape, a ROM or integrated circuit, a magneto-optical disk, a radio or infra-red 

10 transmission channel between the computer module 102 and another device, a computer 
readable card such as a PCMCIA card, and the Internet and Intranets including email 
transmissions and information recorded on websites and the like. The foregoing is merely 
exemplary of relevant computer readable mediums. Other computer readable mediums 
may be practiced without departing from the scope and spirit of the invention. 

15 The methods of Fig. 3A or 3B may alternatively be implemented in dedicated 

hardware such as one or more integrated circuits performing the functions or sub 
functions of Figs 3A or 3B. Such dedicated hardware may include graphic processors, 
digital signal processors, or one or more microprocessors and associated memories. 

The foregoing describes only some embodiments of the present invention, and 

20 modifications and/or changes can be made thereto without departing from the scope and 
spirit of the invention, the embodiment being illustrative and not restrictive. 
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The claims (tefining the invention are as follows: 

1 . A method of demodulating a real two-dimensional pattern, the method 

comprising the steps of: 

5 estimating a quadrature two-dimensional pattern from said real two-dimensional 

pattern using a two-dimensional spiral phase filter; and 

creating a demodulated image by combining said real two-dimensional pattern 
and said estimated quadrature two-dimensional pattern. 

10 2. A method as claimed in claim 1, wherein said estimating step includes 

the steps of: 

generating a frequency domain signal from said real two-dimensional pattern 
using a linear transform; 

applying said two-dimensional spiral phase filter to at least part of said 
15 frequency domain signal to provide a filter signal; 

generating a spatial domain pattern from said filtered signal using an inverse of 
said linear transform; and 

extracting from the spatial domain pattern an estimate of said quadrature two- 
dimensional pattern. 

20 

4 

3. A method as claimed in claim 1, wherein said estimating step includes 

the steps of: 

convolving said real two-dimensional pattern with a complex function to provide 
a convolved spatial domain pattern, said complex function being an inverse Fourier 
25 transform of said two-dimensional spiral phase filter; and 

extracting from the spatial domain pattern an estimate of said imaginary two- 
dimensional pattern. 
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4. A method as claimed in any one of claims 2 to 3, wherein the extracting 

step further includes determining an approximate orientation angle P of at least one fringe 
pattern in said spatial domain pattern. 

5 5. Apparatus for demodulating a real two-dimensional pattern, the 

apparatus comprising: 

means for estimating a quadrature two-dimensional pattern from said real two- 
dimensional pattern using a two-dimensional spiral phase filter; and 

means for creating a demodulated image by combining said real two- 
10 dimensional pattern and said estimated quadrature two-dimensional pattern. 

6. Apparatus as claimed in claim 5, wherein said means for estimating an 
imaginary two-dimensional pattern includes: 

means for generating a frequency domain signal from said real two-dimensional 
15 pattern using a linear transform; 

means for applying said two-dimensional spiral phase filter to at least part of said 
frequency domain signal to provide a filter signal; 

means for generating a spatial domain pattern from said filtered signal using an 
inverse of said linear transform; and 
20 means for extracting from the spatial domain pattern an estimate of said 

quadrature two-dimensional pattern. 

7. Apparatus as claimed in claim 5, wherein said means for estimating an 
imaginary two-dimensional pattern includes: 

25 means for convolving said real two-dimensional pattern with a complex function 

to provide a convolved spatial domain pattern, said complex function being an inverse 
Fourier transform of said two-dimensional spiral phase filter; and 

means for extracting from the spatial domain pattern an estimate of said 
imaginary two-dimensional pattern. 
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8. Apparatus as claimed in claims 6 or 7, wherein the extracting means 
further includes means for determining an approximate orientation angle p of at least one 
fringe pattern in said spatial domain pattern. 

5 

9. Apparatus for calculating a quadrature conjugate of a real two- 
dimensional code pattern, the apparatus comprising: 

a first spatial light modulator for modulating coherent light to produce said real 
two-dimensional pattern; 
10 a first lens for Fourier transformation of said real two-dimensional pattern to 

produce a first Fourier transformed image; 

a second spatial light modulator or spiral phase plate for phase modulating said 
first Fourier transformed image to produce a phase modulated image; 

a second lens for Fourier transformation of said phase modulated image to 
15 produce a second Fourier transformed image; 

a third spatial light modulator or spiral phase plate for phase modulating said 
second Fourier transformed image to produce said conjugate of said real two-dimensional 
pattern; and 

an image sensor for capturing said conjugate of said real two-dimensional 

20 pattern. 

10. A method of demodulation a real two-dimensional pattern, said method 
being substantially as described herein with reference to the accompanying drawings. 

25 11. Apparatus for calculating a conjugate of a real two-dimensional pattern, 

said apparatus being substantially as described herein with reference to Fig. 10 of the 
accompanying drawings. 
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Appendix 1 of specification 
"Spiral Phase Demodulator for Two-dimentional Fringe Patterns" 



m 

Generalised Fringe Pattern Analysis using a Spiral Phase 

Transformation 

Kieran G. Larkin 
Donald J Bone* 
Michael A. Oldfield 

Canon Information Systems Research Australia Pty Ltd. 
*CSIRO Division of Mathematics and Information Sciences 

Abstract 

We introduce an operator with properties that are consistent with the ideas of quadrature and Hilbert 
transformation in two-dimensional signals or images. The 2-D quadrature algorithm presented 
relies on a whole new class of operators we call "spiral phase" or vortex operators, which have 
intriguing orientational properties yet maintain strictly flat spectral response and scale inavariance. 
We demonstrate the remarkable analytic properties of our spiral phase transformer when applied to 
bandlimited images such as fringe patterns and suggest wide applicability for general 2-D Hilbert 
analysis previously considered difficult or even impossible (eg instantaneous 2-D frequency 
estimation, isotropic texture analysis/synthesis, automatic calibration of phase-shifted 
interferograms). The spiral phase transform can be efficiently implemented as a compact 
convolution kernel or as a Fourier space filter. We note that most previous attempts by engineers 
and physicists to extend the Hilbert transform to 2-D have failed to attain isotropic performance. 
Connections between our spiral phase operators and the Cauchy-Riemann equations as well as the 
Riesz transform of harmonic theory are investigated. 
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,1 Introduction 

The concept of analytic (or holomorphic) signal was introduced to communication theory by Gabor 
in 1947 1 for one dimensional (1-D) signals. An analytic signal can be shown to consist of 2 parts; 
the real part is the base signal and the imaginary part is the Hilbert transform of the real part. 
Unfortunately the concept does not naturally extend to 2-D. A number of ad hoc definitions of the 
2-D Hilbert transform have been proposed 2 ~ 8 , all of which have varying degrees of directionality 
(non-isotropism) associated. Typical definitions have quadrant based symmetry. 9 ' 10 In 
multidimensional analysis this extends to orthant symmetry. 1 1 More recently Havilicek et al 12 
have proposed the "analytic image" based on work on 2-D texture analysis and 2-D demodulation 
of AM-FM modelling of image characteristics. A number of improvements are made therein but 
the definition is not isotropic, as illustrated in their test images. Another new development is the 
idea of extending the complex analysis of the Fourier transform to hyper-complex numbers. The 
concept has been called the "quaternionic (W Hamilton 1865) Fourier transform" by Bulow 13 * 14 , 
and it allows an unambiguous definition of the analytic image in 2-D. Unfortunately and 
surprisingly, perhaps, this definition too has a degree of directionality which is apparent in the 
demodulated envelope patterns. 13 Interestingly the same idea to use quaternions and even 
octonions was proposed by Craig in 1996 15 . A number of possible definitions of the 
multidimensional Hilbert transform are considered therein. The quaternionic approach introduces 
an additional phases in the definition of the analytic image; a direct consequence of hyper-complex 
numbers with 3 components. 

Research in geophysics has produced a number of possible definitions of the 2-D Hilbert transform 
over the years. 15-19 The idea of a 2-D signum function necessary for 3-D magnetic field analysis 
was introduced by Nabighian 17 and later applied by Moon et al. 18 However this definition reflects 
the underlying aim to estimate vertical field derivative from two horizontal field derivatives, so 
their definition results in a two component complex vector. Our approach effectively shows how a 
definition of a complex scalar (not vector) 2-D Hilbert transform can be achieved. One of the 
possible quaternionic Hilbert transforms considered by Craig 15 is, in fact, the "classical" case of 
Nabighian's. 

A study of research in the pure mathematical subject of Harmonic analysis, in particular singular 
integrals reveals some interesting parallels. Essentially, the generalisation of the Hilbert transform 
from 1 -dimension to n-dimensions is known as the Riesz transform (named after Marcel Riesz 20 ). 
This may explain some of the confusion in recent years regarding the existence of the Hilbert 
transform for n dimensions, where n>l . In a textbook by Stein 21 the proofs of existence, 
convergence etc of the Riesz transform are derived (see page 56 especially). The Riesz transform 
is, in fact, a Calderon-Zygmund singular operator. 22 Calderon and Zygmund investigated the 
higher dimensional analogues of the conjugate function mapping (or Hilbert transform) in the 
1950's. More recently research into the so called "bilinear Hilbert transform" has been stimulated 
by a conjecture of Calderon's in 1951. However, the bilinear Hilbert transform 23 is actually the 1-D 
Hilbert transform of a function that can be represented as a convolution of two functions, and so 
very different to a Hilbert transform in 2-D. 

Another connection referred to by Nabighian is the generalisation of the Cauchy-Riemann 
conditions to higher dimensions. 24 This connection is also explored by Stein 25 alongside the Riesz 
transforms and the Hilbert transform. Adding to the overall confusion Reich 26 gives a rather 
different definition of a 2-D Hilbert transform based on the work of Calderon and Zygmund 27 . 

We will outline the major properties and derivation of the 2-D spiral phase transform in section 2 so 
that it can be compared with the other transforms such as the Riesz. It should be noted that our 
derivation of the 2-D Hilbert transform analogue is from a very different perspective, and ultimately 
our transform is not quite the same beast! 



we propose a transform which captures the principal properties of the 1-D Hilbert transform and 
maintains a kind of isotropy in 2-D, yet can be defined within the structure of complex functions in 
2-D, without the need for 2-vectors. As we shall see, a direct consequence of this is that these 
operators may be implemented in optical systems (using phase filters) as well as conventional 
complex image processing systems. A simplification of the rigorous transform definition allows an 
easily calculated 2-D Hilbert transform with simple n phase-wrap anistropy but maintains amplitude 
isotropy. 

In addition we define a new class of complex, isotropic, flat spectrum (pure phase) orientation 
operators which we utilise in transform implementation. Our inspiration, in part at least, for these 
operators has been the recent work on spiral phase 28 " 44 (also known as phase singularities, phase 
vortices, spiral singularities, spiral dislocations, and residues). These phenomena occur naturally in 
electromagnetic fields whenever there is a zero in the field intensity and can be characterised by a 
2-D complex Taylors (or Laurent) series expansion about the zero point. 45 It transpires that these 
phase spirals have fascinating propagation properties (refs?) which have been the subject of much 
investigation (refs). In nonlinear media phase singularities can exhibit soliton behaviour 46 . In the 
study of fringe patterns, the common structures of ridge bifurcation and ridge termination can be 
modelled simply by underlying singularities. 47 - 49 In the important area of 2-D phase unwrapping, 
spiral phases are at the heart of ambiguous unwrapping related to noise and other distortions 45 * 50 . 
But our real interest here has been the effect of optical elements with spiral phase structure and their 
effects on EM fields. 29 ' 51-53 In the classic 2 -F optical Fourier transform processor a spiral phase 
can be considered as a complex operator. This turns out to be a very useful way of considering the 
properties of the 2-D Hilbert transform - as a complex multiplier (operator) in the Fourier domain. 
In section ? we will cover some relevant background to the spiral phase operator approach to 
Hilbert transformation. 

The isotropic Hilbert transform opens up new applications in pattern analysis previously considered 
difficult 54 or impossible. 6 * 55 In fringe analysis a number of intractable problems become almost 
trivial because complicated difficulties in spectral zone selection disappear. Many practitioners in 
the area of signal processing 10 are constrained by suboptimal definitions of the multi-dimensional 
Hilbert transform. The spiral phase formalism for the two-dimensional Hilbert transform developed 
in this paper marks a conceptual change from the incumbent linear signum functions to rotary 
signum functions. The proposed change is not unlike the preposterous idea of the screw propeller 
in the nineteenth century, a time when the paddle steamer was dominant. 

2 Definition of the 1-D Hilbert Transform 

According to Ahmed Zayed 56 the Hilbert transform was so named by G.H. Hardy zfrter David 
Hilbert (1862-1943), who was the first to observe the conjugate functions now known as a Hilbert 
transform pair 57 . Major contributions to the theory of Hilbert transforms were made by 
Titchmarsh 58 and Gabor. 1 Following its use in the radar theory of Woodward 59 the applications to 
communication theory developed. 60 " 62 

Perhaps the best place to start is Bracewell's 63 definition of the Hilbert transform in 1-D: 

"{/(*)} =^*/W 

Where * represents the convolution operator. By the Fourier convolution theorem: 



F{«{/W}} = f{^}f(u) 



- 7-) ~~ 



With the forward and inverse Fourier transforms defined as 

F {/(*)} = F (") = $~f(x)exp(-2mtx)dx, 
and 

F ~'{ F ( M )} = /(*) = ^F{u)cxp(+2nux)dx. 
We find 

F{^} = /sgn(«) = 



+i,M> 0 
0,W = 0 
-/,M<0 



In Fourier space the Hilbert transform is obtained by the operation of multiplying by the sign of the 
frequency ordinate u: a very simple operation. Hence frequency space operations are often 
preferred for efficient computation. A crucial property of the Hilbert transform is the spectrally 
independent (that is for any value of the spectral scale factor //) 90° phase shift, or quadrature 
property: 

7/{sin(/#)} = cos(px) 
H{cos(/jx)} = - sin(//x) 

Also f(x)-iH{f(x)} is known as the (complex) analytic signal, and can be represented as 
magnitude and phase: 

f(x)-iH{f(x)} = a{x)e-*« x) = a(x)cosyr(x)-ia(x)smyr(x) 

Again, in Fourier space the situation simply depicted as a function with no negative frequency 
components: 

0,w<0 
F(u),w = 0 
2F(w),w>0 



F{f{x) - iH{f{x)}} = F{u) + sgn( W )F( W ) = 



The above property was used to great effect in the 1980s to revolutionise the analysis of fringe 
patterns, 64 * 65 especially optical interferograms, which were previously analysed manually. The 
initial methods relied on purely one-dimensional analysis, repeated line by line for 2-D patterns. 
However, it was soon discovered 66 that the extension to full 2-D allowed much better 
discrimination between neighbouring frequency components. 
The underlying method for 1-D fringe pattern is as follows: 

f(x) = a(x) + b{x) cos(27rw 0 jc + <p(x)) 

The complex exponential form is: 

f(x) = a(x) + ^^-[exp(2^w 0 x + i$(x)) + exp(-2;nu 0 * - /0(jc))] 
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Courier transtorming the above 
P(u)= A(u) + ^[c(u-u Q ) + C*(u + u 0 )] 

where we define the combined amplitude and phase modulated signal c 

c(x) = b(x)exp(i<p(x)) 

and the corresponding bandpass spectrum 

C(u) = £Jc(x) cxp[-2mux]ix 

It can be shown that, if b(x) and <p(x) are slowly varying functions, then their Fourier transforms 
are spectrally localised. A similar argument applies to a{x). The argument can be made quite 
rigorous, in terms of constraints upon certain derivatives of the functions. Typical functions and 
transforms are shown in Fig 1 below: 




Figure 1(a)- Modulated fringe pattern 
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Figure 1(b). Spectrum of modulated fringe pattern 



Now it is clear that the Fourier transform (FT) has separated the three components of the original 
modulated pattern. The FT can be filtered to remove the central (DC) peak and one of the 



sidelobes leaving just C(u - m 0 ), for example. If we now inverse transform the filtered spectrum we 
get the following: 

frrM ( x ) = C C ( W ~ «o)cxp[+27tfMx)& = cxp[+2mu 0 x]c(x) = exp(27riw 0 x + i<f>(x)) *' 
The next formal step is to find the complex logarithm of the above function: 
log{/™ (*)} = log{^(^)} + i(2nu 0 x + 

In practice the modulus and phase terms are computed as follows: 
nndU^+Kx)) = arctan||^^-| 

We will not consider the further processing of the magnitude and phase information here, but 
instead will consider the basic transformations that have occurred: 

More specifically 

b(x) cos(2mu 0 x + i<p{x)) — > b(x)[cos(2mu 0 x + i<p(x)) + i sin[2mu 0 x + i<p{x))\ 

We see that the process essentially involves the generation of a sinusoidal (quadrature) component 
from a cosine. This is effectively the Hilbert transform of the original fringe pattern, the process is 
obscured slightly in this instance by the removal of the near-DC component of the signal. The 
process could be summarised as a short algorithm: 

1. Remove DC from signal 

2. Fourier transform signal 

3. Multiply by signum(frequency) 

4. Inverse transform 

5. Add this signal (imaginary) to the original function (real and with DC removed^ 
Typically this is the algorithmic starting point for the extension of the Hilbert transform to 2-D, as 
we shall see. 

3 The 2-D Signum function 

Unfortunately the Hilbert transform does not extend naturally into 2-D. Conventional approaches 
typically involve the definition of a 2-D signum function. The simplest approach has been to define 
it as a signum operator of one or other of the frequency coordinates: 
S(u,v) = sgn(w) 
or 

5(w,v) = sgn(v) 
or more generally 

S(u,v) = sgn(w') where u is a rotated frequency coordinate. 

This simple approach, which is inherently one-dimensional, has been used by a number of 
researchers. 2 * 5 « 7 » 8 - l9 » 67 Another common approach has been the alternating quadrant, 3 



S(w,v)= sgn(w)sgn(vj 

Some researchers have proposed a transform based on a number of alternating sectors 4 
S(w,v) = sgn(sin(n0)) 

where n is an integer and <f> is the polar angle in frequency space. This type is a generalisation of 
the preceding definitions. An illustration of the definitions so far is given in Figure 2. 



Figure 2. Signum functions for n=l, 2, 3. 

The main problem with these definitions is the spatial anistropy and the sharp discontinuity between 
neighbouring spatial frequencies. A very good way to test for anistropy is to use a test pattern with 
circular symmetry and apply the candidate Hilbert transform and look for asymmetry in the output. 
From this point of view all the above candidates fail as Figure 3 shows. 













Radially symmetric test pattern 


Halfplane Hilbert transform 


Quadrant Hilbert transform 



Figure 3. Test image outputs for various Hilbert signum function definitions. 



The first pattern in figure 3 is a greyscale representation of a real symmetric function with a 
chirped frequency and an annular amplitude modulation. The two attempts at 2-D Hilbert 
transformation show strong direction artefacts. Extending the number of sectors in the signum 
function merely increases the number of directions, but does not remove the pattern alternation 
between directions, as the various 2-D quadrant functions proposed by numerous researchers over 
the years. 3 - 9 > 1 1 • 68 show. 

The 1-D operator ensures that the functional value ant any frequency is the negative of its value at 
the conjugate frequency (that is to say its 1 80° rotated frequency). The operator is discontinuous at 
a single point, the frequency origin. Ideally we would like something in two dimensions which 
gives a sign change for conjugate frequencies yet avoids discontinuities between neighbouring 
frequencies. Another property we would (probably) also like to maintain is the non preferential 
treatment of different frequencies. We can apply these constraints more formally to an idealised 2- 
D signum function. 

S(m,v) = -S(-«-v) 



and 

\s(u t v) = l| => S(u y v) = exp[iyr(u y v)] 
imply that 

yr(-u -v) = yr(u,v) + 7r 

The constraint upon the neighbouring frequencies is more difficult to state exactly, but it is related 
to the smoothness of the function and the number and type of discontinuities in S. We will limit our 
search to the initial problem space of complex functions of two (real) variables u and v, the spatial 
frequencies. At this point we note that all the recent progress in defining an isotropic Hilbert 
transform have used either complex vector functions of two real variables 17 or quaternions 13 ' 15 » 17 > 
18 to allow some extra degrees of freedom in the problem in pursuit of a solution. It is interesting to 
note that the 2-D Hilbert transform of Nabighian is a complex 2-vector function, and it is not clear 
how this relates to the original complex scalar function, except when applied to potential function 
derivatives (in which case a scalar product with a vector gradient produces a scalar!). 

Figure 4 defines some important quantities in our problem 
space 




Fourier 
Transform 
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Figure 4. Space and spatial frequency coordinates 



Now a possible constraint on S is that is varies uniformly in the radial and tangential spatial 
frequency coordinates q and <p 
dS 

— const 



dr 
and 
dS 



= const 



dip 

Nabighian 17 extends a 2-D potential field problem to 3-D and in the process extends what we term a 
1-D Hilbert transform to 2-D. The problem he solves seems to naturally extend from 1-D to 2-D, 
but only with hindsight! He defines a vector Hilbert transform signum function as follows (in our 
notation): 



\ iu iv •( u . v 1 

y ( M ' v = " / a 2 e -" / 2 2 e y = H~ e - + "^ 

where e x and e y are unit vectors in the x and y spatial frequency directions (u,v) respectively. Our 
approach to the problem has been quite different, and perhaps less direct. We propose a signum 
function as follows: 

S(u, v) = exp(i01 = cos0+/ sin <p — 

Note that the modulus is exactly 1 (except at the origin where it is zero by the Dirichlet condition). 
Our aim in this instance was to produce a function that is a 1-D signum when evaluated along any 
radial slice. The radial slice is further "encoded" with a pure phase function that gradually and 
uniformly varies with the polar angle of the radial slice. The gradual variation is required to avoid 
the sudden sign changes associated with previous 2-D signum definition when evaluated as radial 
functions. Our choice gives the following polar derivatives: 

^ = 0 

dr 

j_as = . 
sd<f>~ 1 

A few observations about our candidate 2-D signum function are in order at this point. Firstly, the 
function is not strictly isotropic, but it does have very strong circular symmetric properties. 
Secondly, the function is what we might call a spiral phase function, or a vortex phase, also known 
as a spiral discontinuity or (spiral) phase singularity, or a complex residue. The numerous names 
reflect the importance of this function in many diverse areas of research. From our present 
perspective the connection between spiral phase and optical diffraction theory, interferometric 
phase-unwrapping and speckle inteferometry is noteworthy, however we do not have the space 
(unfortunately!) to explain all the intricate and fascinating interrelations. 

We are now about half way through the process of defining an isotropic Hilbert transform - we 
have a 2-D signum function that fulfils a number of our criteria. But it is no longer possible to 
define a Hilbert transform purely in terms of real space convolution (frequency space multiplication 
with a signum function) - an extra operation is required to restore isotropy. 
Interestingly, pure mathematicians define an extension of the Hilbert transform to n dimensions 
which they call the Riesz transform. 21 Just two properties are required to define it: The first is 
translation invariance, the second is scale invariance. Our definition actually satisfies these 
constraints because: 1) we have a multiplicative signum function (in the Fourier doifiain), and 2) our 
signum is independent of q. 

4 Demodulation of 2-D Fringe Patterns 

It is useful to consider the stages a test image goes through as we attempt to compute its 2-D Hilbert 
transform. We will avoid the distinction between continuous and discrete data for the time being to 
avoid clouding the issue. Figure 5 shows our previously used test function. We have generated 
an "analytic" image, insofar as we defined a radial symmetric amplitude and a radial symmetric 
phase function 

g(x y y) = b(r)exp[i<f>(r)] = g R +ig g 

We chose a simple, smooth and positive envelope function. The phase consists of a linear function 
of r. By doing this we already have what we believe to be the exact Hilbert transform of the real, 
part in the imaginary part. 





Figure S.Test function for 2-D Hilbert demodulation tests. 



In demodulation problems we can that the cosine component of the signal is known and that the 
quadrature component is to be estimated. In figure 6 various stages in the Hilbert transform process 
are shown. Initially (figure 6(a)) the real input function is Fourier transformed, to produce a real 
symmetric function (by symmetry of FT operation). The Fourier transform shows strong 
localisation in a ring of frequencies and a ± dipole structure is apparent. A recent paper by 
Amidor 69 explores the details of such radially periodic functions and their Fourier spectra. The 
structure if the spiral phase signum function is then shown in terms of real imaginary and phase 
components. The first stage of our Hilbert transformation is the multiplication of the test image 
Fourier transform by the spiral phase function. The result of this is shown in terms of its real and 
imaginary components. In figure 6(b) the product is inverse transformed and displayed as real and 
imaginary or magnitude and phase components. These two components, expressed as a vector 
correspond to Nabighian's definition 17 of a 3-D Hilbert transform (given in his equation 22b). 
Intuitively our expectation of a 2-D Hilbert transform is H{g R } = -g { , or more specifically 



H{b{r) cos^(r)]} = -b{r) sin|>(r)] 



We expect an entirely real function with the same (circular in this case) symmetry as the original 
function. Clearly this is not the case, unless we use the magnitude of the output function (figure 
6(b)) |fc(r)sin[0(r)]|. Using the positive magnitude introduces a feature well known to researchers 

using FTM - phase ambiguity. Figure 7 shows the theoretical analytic image compared to the 
analytic image derived from our spiral phase demodulation followed by a magnitude based 
extraction of the Hilbert component. The phase function looks rather different due to the phase 
ambiguity. There is also a seeming mismatch in the magnitude of the estimated analytic function. 
In fact this is likely to be due to the discontinuity in the radial phase function at the origin, and it is 
not really clear how the analytic function should behave because it can no longer b& considered 
bandlimited. 
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Real test pattern g R {x,y) 



Zero imaginary component 



FFT of test pattern (real) 
G R {u,v) 



Zero imaginary component 





Real part of signum function S(u,v) 



Imaginary part of S(u,v) 




Unit magnitude 



Phase part of S(u,v) 





Real part of S(u,v)G R (u,v) 



Imaginary part of 
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Figure 6(a), Spiral phase signum function operating on test image 
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Imaginary part of 
s(x,y)**g R (x,y) 




Magnitude part of 
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Phase part of 
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Figure 6(b). Output from spiral phase signum operator 
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Figure 7. Analytic Images: theoretical versus spiral demodulated 
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The phase ambiguity appears as a saw tooth profile in the range 0 to n rather than -n to +n as shown 
in figures 8(a) and 8(b). Mathematically it arises from the arctangent function returning the 
principal value in the correct quadrant of the full 2n range. 
sin(^) j_ 



arctan 



co: 



<+). 



mod 2/r [0 + 7v]-7V 



arctan 



|sin(^)| 1 = f mod^], if mod 2jr [*] < n 
cos(0) J \k - mod n [<f>\ if mo& 2n [<fi\ > Tt 




Figure 8(a). Modulo 2n phase calculation 
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Figure 8(b). Phase calculation with n ambiguity 



We shall show later how the phase ambiguity problem may be overcome for many cases of interest. 
5 Compensating for the spiral phase encoding 

We saw in the previous section how we could define a smoothly varying Fourier space signum 
function with a spiral phase structure. The simple output from this filter consists of two separate 
components which do not immediately fit in with our preconceptions of what a 2-D Hilbert 
transformer should do. To understand more clearly it is helpful to consider a pattern with just one 
pairs of conjugate frequency components in the Fourier domain. Figure 9(a) shows localised fringe 
pattern structure that gives rise to significant energy in the sidelobes shown in Figure 9(b). 



Local fringe pattern 




Figure 9(a)* Local fringe pattern structure 

In Figure 9(a) the 2-D fringe pattern located around coordinates (x 0 ,yo) with spatial frequency 
parameters (u,v) and amplitude b is 

f{x,y) = a + b co^2tv{u 0 {x - x 0 } + v 0 {y - y 0 }) + v] 

Here we assume that the frequency components are slowly varying so that they are effectively 
constant over the small region of interest, Q , around (x 0 ,yo). The fringe offset and modulation, 
a(x,y) and b(x 9 y) are also assumed to be slowly varying functions. 
The Fourier transform (FT) of the above fringe pattern region is 
F n( w > v ) = JJ / (*» y) exp[-2fltf(ux + vy)]dxdy 

= A n (w,v) + 5 n (w-w 0 ,v-v 0 ) cxp[ / yr] + B a ( u + u 0 , v + v 0 ) exp[-/ y/\ 

Figure 9(b) shows the sidelobes with the central (near DC) lobe suppressed. 




Spectral lobes 
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Figure 9(b). Spectral lobes from localised fringe structure 

The local spatial frequency of the pattern satisfies the following equations: 

«o+ v o =°l 



Ik 



u 0 

The normal to the fringes being at an angle f3 Q to the x-axis. The conventional Fourier transform 
method (FTM) of fringe analysis is based on the assumption that the angle (3 Q is confined to a small 
range (certainly much less than n radians). More generally fringe patterns with fringes have the 
orientation angle /? 0 covering the entire range from 0 to 27r radians. 

If the 2-D signum function (phase spiral) is applied to the pattern spectrum then we get a sign 
reversal in the lobes: 

G( u , v) = [F n ( u , v) - A( u , v)] exp[i y/] = B a ( u - u 0 , v - v 0 ) exp[+i>] exp[+^] 

+ B n (u + u 0 , v + v 0 ) cxp[-iy/] cxp[-/ <p] 

Now we can transform back (implicitly assuming that <p{u,v) « ±f3 0 remains constant over each 
lobe). This is equivalent to assuming that the lobe subtends a small polar angle. 



g{x, y) = £TJ^ G ("» v)exp[2/n(wjc + vy)]dudv 

= -i exp[i/3 0 ]b s'm[lm(u 0 {x -x 0 } + v 0 {y - y 0 }) + y] 



So all that remains to be done to recover the sine component of the local fringe pattern is to 
multiply by exp[— 1/? 0 ], The process is repeated for all fringe regions with their corresponding 
fringe orientation angles p. 



exp[-/#)]£(x, y) = -ib sin[2^*(w 0 {x - x 0 } + v 0 {y - y 0 }) + y/] 



We conjecture that if fringe pattern has suitable constraints upon the smoothness of the fringe 
parameters (orientation, spacing, offset, and modulation), then the above equation is a good 
approximation and gives the quadrature component to the original cosinusoidal function. 
It was actually shown earlier that the phase function is implicit in the function g(x,y) 9 but with a k 
phase ambiguity due to the loss of sign information in a square root operation. One way to view 
this is in terms of the fringes in figure 9(a). If the fringes were rotated 180°, then the diagram 
would be completely unchanged. However, there are certain constraints upon how fast, and in what 
way, fringe orientation will vary in many instances. In which case it may be possible to use these 
constraints to resolve the n ambiguity. As we shall show, it is easier to apply such constraints to the 
orientation map than it is to apply to the final fringe quadrature estimate. Typically the orientation 
map fi[x,y) can be estimated by considering the gradient of the fringe map and extracting the 
fringe normal. This is similar to the work of Nabighian who computes the z derivative of the field 
from a scalar product of a 2-D Hilbert transform vector and a gradient vector. In our case this 
approach can give the same ambiguous results as shown in figure 7 (if we take into account that we 
are not interested in derivatives of our function, but the functions themselves). It turns out that we 
can easily overcome the ambiguity by varying degrees by controlling our estimate of y) . 
Firstly it is important to have estimates of P{x,y) even when |Y/(x,y)| = 0, which occurs regularly 
on the peaks and troughs of the fringes. If these estimates are not available we get very poor 
results. So we suggest using a method based on 2-D energy operators 70 which include second 
derivatives and give uniformly distributed estimates of energy and orientation across modulated 
sinusoids. Other estimates are also possible, but such an investigation would be a major work in 
itself. We assume a reliable orientation estimator is available. The energy based "gradient" operator 
we actually use gives estimates of 2/?. Using a convenient complex notation we get an estimate of 
exp[2//?(;c,;y)] to be precise. Such complex notation is very useful in that it inherently wraps the 
estimate exactly in line with our inherent (n) ambiguity. In particular, the estimate of /3 involves 
first finding the square root of the complex exponential then the argument. 



If we apply the implied algorithm and just take the positive root for simplicity we get a rather nice 
result for simple closed patterns as is shown in figure 10. Most of the ambiguity induced phase 
discontinuities are now gone and we have a modulo 2n phase with an overall n "flip" vertically 
down the centre. For comparison the conventional FTM 1-D Hilbert approach is also shown. 




(b) Estimated analytic image using 
simple orientation map (mag/phase) 
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(c) Estimated analytic image 
conventional halfplane Hilbert 




(d) Estimated analytic image using 
unwrapped orientation map 
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Figure 10. Analytic Images: Orientation map based demodulation. 
6 Compensating for ambiguities in the orientation estimator 

The estimate of the 2-D Hilbert transformed image in the previous section still has one unwanted 
feature which is related to an overall ambiguity in the fringe direction resulting in ^loss of sign 
information for our Hilbert transform. There is an elegant solution to this problem which is related 
to the physical constraints upon the expected fringe pattern structure. If we include a phase 
unwrapping process with the inherent square root operation on the complex orientation map we can 
avoid sharp discontinuities in the (3 estimate due to inappropriate choice of the sign. Essentially we 
apply an a priori restriction on the smoothness of /?. There are many ways to do this, but we are 
fortunate in our approach to the problem because of the strong topological constraints on orientation 
maps of fringe patterns. Previous work on the analysis of oriented patterns 47 * 49 ' 7 1 used the 
orientation of the patterns to help extract useful pattern descriptors. 

A number of recent works on phase unwrapping have summarised the performance of panoply of 
algorithm developed in the last few decades. In particular the special processing of phase 
singularities in phase maps using complex variable analysis (residues etc) has matured. These 
methods can be used to advantage to unwrap our pure phase orientation map. In our particular case 
we can have singularities of order +2, +1,-1, whereas conventional interferograms (eg Synthetic 
Aperture Radar interferograms, speckle interferograms etc) only have +1, and -1, with higher 
orders being so unlikely as to be virtually non-existent. 



We can view our overall process as being optimised to concentrate the unwrapping step (which is 

needed in-ALL existing fringe pattern- analysis techniques) in the most easily solved part of the 

process: the orientation unwrapping. This was not a prior requirement but a natural outcome of our 
process. The work of Nabighian only really considers functions whose DERIVATIVES obey the 
Hilbert transform relation. The functions themselves have VECTOR Hilbert transforms in his 
analysis. We have essentially taken the extra step to reduce the 2 vector Hilbert transform to a real 
function, in line with our pre-conceptions of what a 2-D Hilbert transform should be. 

7 Singular Integrals and the Riesz Transform versus the Ouaternionic Hilbert 

Transform 

Subsequent to our discovery and implementation of the spiral phase (vortex) "2-D Hilbert" 
transformation, our literature searches led to some work on singular integrals by Calderon and 
Zygmund. It transpires that the second order Riesz transform is analogous to the proposed 2-D 
Hilbert transform of Nabighian 17 and the first quaternionic transform example by Craig. 15 Stein 21 
gives the following example of a n-dimensional Riesz transform in terms of its (n) convolution 
kernels: 




j = 1 — » n 



In terms of the multipliers (viz the Fourier transforms of the kernels): 
m A u ) = ^ j=l~>n 

So the multiplier effectively defines an n-dimensional signum function, although it is not typically 
viewed in this way by the pure mathematicians. An important point to note again is that such a 
transform maps a 2-D real scalar function to a 2-D real vector (2 component) function, whereas are 
proposal is to map a 2-D real scalar function to a 2-D complex function, then to a 2-D real function 
(using orientation information). We feel that such an approach is more consistent with our target 
applications in fringe analysis, image processing and optical diffraction theory, notwithstanding the 
inevitable difficulties of then extending beyond 2 dimensions. Indeed all proposed operations are 
possible within the domain of complex operators acting upon complex images. 
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